function [g_t,delta_g_t,g_p,g_norm,res,pu,delta_pu]=gradient(r,U,delta_U,n)
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

c1=1/c;
a=a_t';
res=zeros(1,n_t);
pu=zeros(1,n_t);
delta_pu=zeros(1,n_t);
for i=1:n
  pu=pu+r(i)*U(:,i)';     %pu=p(u) the peak selection
  delta_pu=delta_pu+r(i)*delta_U(:,i)';
end

[fu,dfu,sfu]=nonlinear_terms(pu);
                                                                                
res=c1*delta_pu-a_t'.*pu+fu;
g_p=assempde(b,p,e,t,c1,a,res);        %The negative gradient of J at p.
g_t=pdeintrp(p,t,g_p);                %The negative gradient of J at t.
g_norm=res*(g_t.*ar)';
g_norm=g_norm^0.5;                    %Its norm
delta_g_t=c*(a_t'.*g_t-res);
